Surface tension of the isotropic-nematic interface 
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We present the first calculations of the pressure tensor profile in the vicinity of the planar interface 
between isotropic liquid and nematic liquid crystal, using Onsager's density functional theory and 
computer simulation. When the liquid crystal director is aligned parallel to the interface, the 
situation of lowest free energy, there is a large tension on the nematic side of the interface and a 
small compressive region on the isotropic side. By contrast, for perpendicular alignment, the tension 
is on the isotropic side. There is excellent agreement between theory and simulation both in the 
forms of the pressure tensor profiles, and the values of the surface tension. 
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The nematic liquid crystal phase is distinguished from 
the isotropic liquid by the existence of long-ranged ori- 
entational ordering of the constituent molecules jl| . The 
structure of the planar interface between coexisting ne- 
matic (N) and isotropic (I) phases is of fundamental in- 
terest: along with the liquid-vapour surface, it is one 
of the simplest fluid-fluid interfaces. Molecular-scale in- 
teractions, capillary-wave fluctuations, and the coupling 
between translational and rotational degrees of freedom, 
may all influence the interface properties. The surface 
excess free energy, or interfacial tension 7, is a key pa- 
rameter. The variation of this quantity with the direction 
of molecular ordering in the liquid crystal, the so-called 
director, dictates the molecular alignment which will be 
preferred by the surface: the phenomenon of anchoring, 
an essential feature of liquid crystal physics j|. Deter- 
mining the surface tension of the N-I interface for a given 
molecular model by computer simulation, and predicting 
it from first principles by theoretical methods, is clearly 
desirable. For a planar interface, 7 may be calculated by 
integrating the difference between normal Pn arid trans- 
verse Pt pressure tensor components with respect to po- 
sition |,| : 7 = dz [P N (z) - P T (z)}. 

In this paper, we present the first discussion of the 
pressure tensor profile through the nematic-isotropic in- 
terface in the framework of a simple density functional 
theory in the Onsager approximation (^|, modelling the 
molecules as freely-rotating and translating rod-like rigid 
bodies. We shall show that there is a significant differ- 
ence between the two situations of most interest, when 
the director is oriented respectively parallel to, and nor- 
mal to, the interface. Previous studies calculat- 
ing the surface tension from free energy differences, have 
indicated that, in the limit of infinite rod elongation, the 
minimum surface free energy is obtained when the rods 
lie parallel to the N-I interface; in other words, this in- 



terface favours planar anchoring. We provide additional 
insight into this result, through an examination of the 
normal and transverse pressure. 

We back up these predictions by carrying out the first 
computer simulations to give the pressure tensor profiles 
for equilibrium N-I interfaces, confined between paral- 
lel walls under different anchoring conditions, using the 
same molecular models as in the theory, Our results 
are further supported by the first computer simulations 
of free standing nematic films in coexistence with the 
isotropic phase. From the simulation viewpoint, this 
study is much more challenging than determining the sur- 
face tension of the liquid-vapour interface fl0|| since the 
densities and other properties of the two phases are very 
similar. To our knowledge, the only other such studies 
have involved two phases maintained out of thermody- 
namic equilibrium 

We employ two very similar molecular models in this 
work, in both simulation and theory. Both models are ax- 
ially symmetric, with molecular orientations represented 
by unit vectors Uj and centre-of-mass positions by vec- 
tors Ti. In the first model, the molecules are taken to be 
hard ellipsoids of revolution of axial ratio k = A/B = 15 
where A is the length of the symmetry axis and B that 
of the transverse axes; the bulk phase diagram of this 
model is well known | fl3[ , and the interfacial structure has 
been recently reported by one of us |3. For hard par- 
ticle models, the temperature enters only trivially into 
the thermodynamics, and we define units of energy by 
setting k^T = 1 in this work. Length units are fixed by 
taking B = 1, mass units by setting particle mass M = 1, 
and the moment of inertia was fixed at I /MB 2 = 50. In 
the second 'soft ellipsoid' model, a continuous pair po- 
tential v 12 = 4e(s 12 12 - s^) + 1, for s 12 < 2 1 / 6 , v 12 = 
otherwise, was used. Here, s\ 2 — {r\ 2 — o\ 2 + B)/B 
is a scaled and shifted separation where the 'diameter' 



1 



012 = f(r 12, Ui, U2; A, B) is a standard approximation 
p5| , p^| to the contact distance of two ellipsoids (again 
with A/B = 15 here) at given relative orientations; 
r i2 = ri - r 2 , r 12 = |r 12 | and f 12 = r 12 /r 12 . The po- 
tential is purely repulsive and varies quite strongly over 
a short range of interparticle separation. We choose an 
energy parameter e = 1, again take k B T = 1, and set 
other parameters as for the hard ellipsoids. 

In planar geometry, with translational invariance in 
the x and y directions, we adopt the Irving-Kir kwood 
convention jlTj so as to write components of the pressure 
tensor in the form fl|It 



P aa (z) = k B Tp(z) 



(1) 



- / dr i2 / dui / du 2 r° 
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x / d( g {2) (r 12 ,z - (z 12 ,z + (1 - <>i2;ui,u 2 ) 
Jo 

a = x,y,z. 

Here, p(z) = Jdug^(z,u) is the number den- 
sity, ^'(zjii) the (position- and orientation-dependent) 
single-particle density, and g^ 2 \r 12l z\, z 2l m, u 2 ) the 
pair density. The normal component P^[z) = P zz (z) — 
P is independent of z throughout the system, while 
the transverse component Pt(z) = P X x(z) = P yy (z) is 
not. Far from the interface, the pressure tensor becomes 
isotropic, Pn(z) = Pt(z) = P. 

For simulations of the soft-ellipsoid system, the pres- 
sure tensor was routinely computed from Eqn (|l]) via 
the intermolecular forces ||,|l0|,|l8| . For the hard-ellipsoid 
system, it was calculated by identifying nearly-contacting 
pairs of ellipsoids Jl^,0 . Such pairs are defined as those 
which would overlap if their dimensions are uniformly 
scaled by a factor A, very slightly larger than unity. Tests 
were carried out with several choices of A: all of the re- 
sults here are obtained with A = 1.01, the smallest value 
to give acceptable statistics, which was also found to be 
satisfactory in ref fiS| ], 

The construction and minimization of the Onsager 
free-energy functional followed closely the procedure de- 
scribed before . The excess part of the Helmholtz 
free energy T = J- ld + J 7ex is expressed in terms of 
«(r,u): 

^ ox [p (1) ] = -~fcsT J dn J d Ul J dr 2 J du 2 (2) 

X£ (1) (ri,ui)0 (1) (r 2 ,u 2 ) /(ri2,ui,u 2 ) . 

where /(ri 2 ,ui,u 2 ) = exp(— vi 2 /k B T) — 1; this is the 
leading term in a low-density virial expansion. Fluctua- 
tion effects, which come in through higher-order terms, 
were neglected. The most important fluctuations are 
capillary waves, which only slightly broaden the profiles 
in small systems, and do not affect the interfacial ten- 
sion. Since g^ is independent of x and y, f may be 



integrated over these coordinates. The logarithm of the 
density is expanded in spherical harmonics, on a dis- 
crete grid in z with an interval 5z — B/5 = A/75: 
lngW^u) = ^i m Ci ; / m y im (u). Periodic boundary con- 
ditions apply in the z direction. In the current appli- 
cation, this expansion was taken to fourth order in the 
Ylm (u) (restricting to second order was found to change 
the coexistence pressure and densities by ~1%). All rel- 
evant orientational integrals were expressed in terms of 
spherical harmonic coefficients of the density up to order 
10 (reducing to order 8 was found to change the coex- 
istence pressure and densities by ~0.01%). The param- 
eters Ci-i m were varied to minimize the grand potential 
= PF[qW] - [3pN of an I-N film system at the 
coexistence value of the chemical potential //, determined 
by a preliminary bulk calculation. 

The normal pressure is constant through the sys- 
tem, and equal to —ft/V in both bulk phases at co- 
existence. The transverse pressure tensor profile was 
calculated directly from the minimized density profiles. 
The definition consistent with Onsager's approxima- 
tion and the use of the Irving-Kirkwood convention in 
the simulations is easily derived from eqn ([!]), by in- 
serting the low-density form g( 2 \ri 2 , z%, z 2 \ Ui, U2) ~ 
g^\zi, Ui)^' 1 ) (z2, U2) exp(— vi 2 /k B T) and carrying out 
an integration by parts. The result is 

P T (z) = k B Tp(z) (3) 

-\ k vT j dr 12 j dUl / dU2 /( r 12> U l' U 2) 

X / dC g (1) (^-^12,U 1 )^ 1 )(z+(l-C)^12,U 2 ). 

Jo 

All simulations were carried out at constant volume. 
A cuboidal simulation box was employed in all cases, 
with periodic boundaries in the x and y directions. Con- 
ventional, microcanonical ensemble molecular dynamics 
(MD) simulations of soft ellipsoids employed periodic 
boundaries also in the z direction, and used a reduced 
timestep St — 0.002. Several systems of N = 7200 parti- 
cles in boxes of size L x /A = L y /A — 2.5, L z /A « 20 were 
studied. At average number densities 0.24 < pAB 2 < 
0.255, from an initially isotropic configuration, one or 
several nematic films were observed to assemble sponta- 
neously, with planar interfaces normal to the z-direction 
and, in every case, the director lying parallel to the inter- 
face. In addition, from an initially nematic configuration 
aligned in the z direction, a similar two-phase system 
nucleated with the nematic director turned through 90 
degrees to lie parallel to the interface, within roughly 1 
million time steps. 

One of the systems at density pAB 2 = 0.255, contain- 
ing a nematic film of width ~ 13-B, was further equili- 
brated over 7 million time steps, monitoring the pressure 
tensor profiles; run averages were then accumulated over 
4.7 million steps. In the absence of walls, the interfacial 
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tension can be calculated simply from the difference be- 
tween the system averages of Pn and Pr ■ This has been 
done independently using the transverse component P yy , 
which happens to be roughly parallel to the director, and 
the transverse component P xx , which is perpendicular, 
giving respectively 1 B 2 /k B T = 0.022±0.01, 0.015±0.01. 
Thus, the free film simulations yielded clear evidence that 
the anchoring at the interface is planar, provided a rough 
estimate of the interfacial tension, and indicated an un- 
expected difference between the transverse stress compo- 
nents P xx and Pyy . To provide more convenient control of 
the director orientation relative to the interface, a more 
systematic study was therefore conducted in a confined 
geometry. 

Both hard and soft ellipsoids were simulated using con- 
ventional constant-JVVT Monte Carlo (MC) j2§ using 
N « 8000 particles in a box of dimensions L x /A = 
Ly/A = 2.22, L z /A = 22.73. Flat walls with a short- 
range orienting field at the box ends allowed nematic 
films of width approximately 8A to be stabilized in a de- 
sired director orientation, separated by a central isotropic 
region of width approximately 7 A. Full details of the 
system, together with a description of the film structure, 
have been given elsewhere In this paper we focus 

on the cases of normal (along z) and in-plane (along y) 
orientation of the director relative to the interface. For 
each model, and each director alignment, we report the 
results of averaging over eight independent runs. For 
the hard ellipsoids, more than 10 6 Monte Carlo sweeps 
were allowed for equilibration, and the averages were ac- 
cumulated over a further 7.5 x 10 5 sweeps (one sweep 
is one attempted move per particle). For the soft ellip- 
soids, the corresponding figures are 3.5 x 10 5 sweeps and 
6.25 x 10 5 sweeps. Profiles of number density p(z), ori- 
entational order parameter S(z) (the largest eigenvalue 
of the second rank order tensor [ H ) and the pressure 
tensor were accumulated. To facilitate the precise com- 
parison between simulation profiles and theoretical pre- 
dictions, a correction was made for the fluctuations in 
the interface position, as described elsewhere Jl4|]. Er- 
ror bars were estimated by assuming a normal distribu- 
tion of the results of the eight separate runs. For the 
in-plane cases of director alignment along y we detected 
(as for the free films) an extremely small discrepancy 
P zz = P xx ^ Pyy in the nematic phase far from the inter- 
face. We ascribe this to correlations of aligned molecules 
across the periodic box: to eliminate this finite-size effect 
from the very sensitive surface tension calculation we set 
Pt(z) = P xx {z) in these cases. For normal alignment, 
we found the x and y directions to be equivalent and set 

Pt(z) = (P xx (z)+Pvv(z))/2- 

Coexistence parameters from theory and simulation 
are compared in Table |. As is well known for 
elongation A/B = 15 the Onsager theory overestimates 
the coexistence pressures and densities by up to 30%. 
Profiles in the vicinity of the interface are compared in 



Figures |l| and g. To help locate the interface, we show 
the number density curves for both simulation and the- 
ory, scaled by the bulk coexistence values pi,pn thus: 
p*(z) = (p(z) — pi)/(pn — pi), and similarly for the ne- 
matic order parameter S*(z). The zero of z is taken at 
the point of inflection of S*(z); as discussed elsewhere 
the density profile lies at z « —A/3 relative to this 
(on the nematic side). In the figures we show both the 
pressure tensor difference Pn — Pt{z) and the running 
integral J z dz Pn — Pt( z ) taken from a point zq in the 
bulk nematic phase far from the interface. The results, 
for the two different director alignments, are dramati- 
cally different. For parallel, in-plane, alignment, there is 
a large tension (low transverse pressure) on the nematic 
side of the interface, followed by a small compressive re- 
gion (high transverse pressure) on the isotropic side, as 
seen in simulations of the liquid- vapour interface of a sim- 
ple fluid p3UTo|| . The surface of tension defined by the 
first moment of Pn — Pr(z) 0] lies on the denser side 
of the interface. For normal alignment, there is a com- 
pressive region on the nematic side, and a region under 
tension on the isotropic side: the surface of tension lies 
on the side of the less dense phase. Such behaviour is 
seen in the idealized penetrable sphere model of the liq- 
uid vapour interface Q but to our knowledge this is the 
first time it has been observed for a molecular model with 
realistic interactions. The step change in the running in- 
tegrals gives estimates of the surface tension, reported in 
Table ffl. We find that the planar orientation is favoured, 
as it has the lower excess free energy; the same result 
was obtained by Onsager theory in the infinite elongation 
limit P|, although the corresponding surface tensions in 
that case (appropriately scaled by A and B) are about 
40% larger than here. 

Although the simulation results are subject to signifi- 
cant error bars, we have succeeded in resolving the vari- 
ation of transverse pressure through the interface, and 
the profiles are very similar, both in form and magni- 
tude, to those predicted by theory. Indeed, the level 
of agreement between theory and simulation results is 
remarkable, bearing in mind the absence of adjustable 
parameters in the Onsager theory. In principle, the loca- 
tion of the surface of tension depends on the convention 
adopted for the transverse stress, e.g. Irving-Kirkwood 
vs Harasima [Q, but in the current study we find that 
the profiles are only slightly dependent on this choice; 
moreover, the results for such molecular shapes should 
be indicative of the behaviour for more elongated parti- 
cles, so we expect them to be rather generally applicable. 
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FIG. 1. Profiles near the N-I interface for hard ellipsoids. Top: surface tension integral J* dz P N — Pr(z). Middle: pressure 
tensor difference Pn — Pt(z). Bottom: reduced density p*(z) and order parameter S*(z). Simulation results are points with 
error bars, predictions of Onsager theory represented as lines. 
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FIG. 2. Profiles near the N-I interface for soft ellipsoids. Notation as for Figure 1. 
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